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A simple optimization scheme is used to compute the density-density response function of an 
electron liquid. Higher order terms in the perturbation expansion beyond the random phase ap¬ 
proximation are summed approximately by enforcing the constraint that the spin density pair cor¬ 
relation functions be positive. The theory is applied to the 3-D homogeneous electron gas at zero 
temperature. Quantitative comparison is made with previous theory and data from quantum Monte 
Carlo simulation. When thermodynamic consistency is enforced on the compressibility, agreement 
with the available simulation data is very good for the entire paramagnetic region, from weakly to 
strongly correlated densities. In this case, the accuracy of the theory is comparable to or better 
than the best of previous theory, including the full GW approximation. In addition, it is found 
that the spin susceptibility diverges at a lower density ( r s « 107) than the current estimate for the 
liquid-solid transition. Application of the theory to inhomogeneous electron liquids is discussed. 

PACS numbers: 71.10.Ca,71.10.-w,71.10.Hf 


I. INTRODUCTION 

In many practical calculations of electronic structure, 
such as for semiconductors and molecules, the density- 
density response function \ plays a central role. The 
popular Kohn-Sham electronic density functional theory 
(DFT)[Mj, and Hedin’s “GW” approximation (GWA) 
of many-body perturbation theory]?] (5] are two such 
methods that use \- 

In DFT, x is an important ingredient in the exchange- 
correlation functional, E xc , which is the main object 
for approximations in the theory. [2] The workhorse lo¬ 
cal density approximation (LDA) and its generalizations 
avoid explicitly estimating \ for the inhomogeneous liq¬ 
uid. Instead, they approximate E xc by relying on know¬ 
ing the correlations of a simpler system, the homogeneous 
electron gas, i.e, jellium. [B] The correlations of this latter 
system have been computed accurately using theory and 
quantum Monte Carlo (QMC) simulation. [7] While be¬ 
ing very accurate for the electron structure of molecules 
and electron density of many systems, DFT-LDA and its 
variants have the limitation of not being able to describe 
accurately band gaps, nor London dispersion (i.e., van 
der Waals) interactions very well. [2] The former is crucial 
to understanding the properties of semiconductors, while 
the latter is important to understanding weakly bonded 
systems, such as water. A recent branch of DFT uses 
the adiabatic-connection fluctuation-dissipation form for 
the correlation energy, which is the troublesome part of 
E xc . In this branch there is much current interest in the 
use of the random phase approximation (RPA) for y.[8| 
Unlike the LDA, this DFT-RPA does predict reasonable 
band gaps for many solids, [9] and, with some extra ef¬ 
fort, reasonable London dispersion interactions between 
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In GWA, x enters through the screened potential W 
to compute ultimately the one-particle Green’s function, 
G, from which the band structure is extracted. HE] In 
perturbative calculations around a DFT reference state, 
known as “one-shot GW,” x is almost always taken to 
have a simple RPA form. In the full version of GWA, 
X is computed in a self-consistent loop along with G. 
Being able to predict accurate band structure for many 
solids, GWA has rapidly become one of the main meth¬ 
ods for computing the electronic structure of crystalline 
materials. [HIS] A drawback to this approach though 
is that the full, self-consistent solution often does not 
improve the predictions of the theory, and can make it 
worse. m Application of the GWA to jellium has shown 
that the main weakness is the expression for x obtained 
from G.|B] 

For both approaches then, a better way to compute x 
would be very desired. In this paper, a simple scheme 
called range optimization is described that goes beyond 
the RPA for x to accomplish this task. Range optimiza¬ 
tion was originally used to improve the RPA theory of 
classical molecules. This “range-optimized” RPA (RO- 
RPA) theory was able to describe very well the equilib¬ 
rium structure and thermodynamics of strongly charged 
polyelectrolyte solutions. used It has also been applied 
successfully to the study of the hydrated electron. p~8] 
Given that a number of weaknesses of the GWA be¬ 
came apparent only after the full version had been im¬ 
plemented for jellium, it is important for the theory here 
to be analyzed first for that system. This is the intent of 
the present paper. As will be shown below, the scheme 
applied to jellium greatly improves the predictions of the 
theory over the basic RPA. A new algorithm to imple¬ 
ment range optimization, valid for inhomogeneous liq¬ 
uids, is also described here. 
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II. THEORY 


Let the analysis be confined to a homogeneous elec¬ 
tron gas, i.e., jellium, in 3-D at zero temperature in the 
paramagnetic phase. As a reminder, jellium has a bal¬ 
ancing background of positive charge that is uniform and 
rigid. [6l IT9l 

The theory here is centered around the time-ordered, 
spin density-density response function (or spin polariza¬ 
tion propagator), Xij( x i x ')i where i and j label the spin 
(t and j,), and x = {r, t}, with r being the position and 
t the time. It is defined as: 

^ J_ (y 0 \T{5hi(x)5h : >(x')}\'I/o} 

in (H'ol'I'o) ■ 1 j 

Here, |d> 0 ) is the system ground state, and T denotes the 
time ordered product. Also, 



Sfii(x) = hi(x) - (hi(x )), (II.2) 


where n, ( x ) is the number density operator in the Heisen¬ 
berg representation for spin i electrons at point x, na and 
( hi(x )) denotes its ground state average. For the param¬ 
agnetic phase, (h-f(x)} = {fi±(x)) = n/ 2, where n is the 
average electron density. 

Jellium is translationally invariant and \ is symmetric 
with respect to x x ', so Xij{ x , x ') —> Xij( r T T ) > where 
r = |r — r'| and t = t — t'. For this case, it is helpful to 
work with the dual Fourier transform, Xijik,uj), k = |k| 
being the wavevector and u> the frequency. 

A number of useful quantities can be obtained from 
Xij(k, w).[5] First is the system compressibility K: 


Ko 

K 


-N( 0) lim 

fc->0 


X(k, 0) 


■ v(k) 


(II.3) 


where y(fc, uj) = JA • Xij(k, uj) is the total density-density 
response function, and N( 0) = mkp/i 7r 2 ?i 2 ), with m be¬ 
ing the electron mass, and = (37i^n) 1 / 3 being the 
Fermi wavevector. [6] Also, Kq is the compressibility of 
the non-interacting gas, and v(k) is the Fourier transform 
of the electron-electron Coulomb potential, v(r) = e 2 /r, 
with e being the electron charge. Second is the spin sus¬ 
ceptibility xs : 


Xp 

XS 


-N(0) lim 

*->• o 


2(XTt(M)-Xt-|.(M)), 


, (II.4) 


where Xp is the spin susceptibility of the non-interacting 
gas. Third are the partial static structure factors: 

c\- j- r OO 

Sij{k) = — / duJXij(k,u), (II.5) 

Ttn Jo 

which are real, and have exploited that Xij(k,uj) is sym¬ 
metric about uj = 0. Fourth are the spin-spin pair corre¬ 
lation (or radial distribution) functions: 


9 ij(r) = 1 + ^ j [Sij{k ) - Sij] exp{ik • r), (II.6) 


where Sij is the Kronecker delta. The pair correlation 
function gij (r) is proportional to the equilibrium proba¬ 
bility density of there being an electron of spin j a dis¬ 
tance r from one with spin i at the same time. As a 
density, g%j{r) is strictly positive for all r. Last, the cor¬ 
relation energy E corr can be obtained from .S',, i /.• i.! 19. 
Two expressions for E corr were used and they are given 
in Sec. nm below. 

Define the matrix inverse of Xij{k,uj) by 
Yj S Xis(k,u))x^(k, u) = Sij. This inverse can be 
represented exactly as: 

Xy{k,w) = X^(k,u) - Vij(k), (IL7) 


where Vjj{k) = v(k ), and Xij(k,uj) is the proper spin 
density-density response function (or proper spin po¬ 
larization propagator). [U QjJl It is helpful to rewrite 
Eq.(|IL7|) as: 


Xi/(*, w) = X 0 ij{k, uj) - Vij(k) - Uij{k, w), (II.8) 

where mj (k, uj) = x^- (k, uj) - X^ 1 (k, uj) and Xoij (k, uj) = 
SijXo(k,uj)/2, with xo(k,uj) being the total density- 
density response function of the non-interacting gas, i.e., 
the Lindhard function. This Lindhard function can be 
computed analytically. m 

Setting Uij(k,w) to zero reduces Eq. (III.8 ) to the famil- 

dLEk 


iar RPA expression for Xij -0 As is well known, the RPA 
tends to work well if the interactions are weak. However, 
for strongly interacting systems it works less well, caus¬ 
ing, for example, the pair correlation functions gij{r) to 
be (very) negative at small r. 

More recent research on jellium has steadily improved 
upon the RPA. Almost all of these theories have worked 
with the one-component expression for the total density- 
density response function, x{k,ui), by developing accu¬ 
rate approximations to the static and dynamic local field 
factors, G+(k) and G+(k,u), respectively. These are de¬ 
fined by: 


X {k,uj) = x o {k,u)-v(k)[l-G+}, 


(II.9) 


where G + denotes G + (k) or G+(k,uj). In this manner, 
the one-component optimized potential, u = —vG +. 

The best of these theories now agree with QMC sim¬ 
ulation data for the paramagnetic state for the corre¬ 
lation energy E corr within a few percent for the den¬ 
sity range of most metals, 2 < r s < 6.[7j [^Dj Here, 
r s = pq/o-o = l/(a/cFao)j where xq = (3/(47m)) 1 / 3 is the 
average distance between electrons, ao is the Bohr ra¬ 
dius, and a = (4/(97r)) 1 / 3 . A limitation of these theories 
though is that they usually apply only to the paramag¬ 
netic, i.e., zero polarization state. This local field factor 
approach can be extended to examine partially polarized 
states, and thus give information about the jellium phase 
diagram. m However, the cost is an increase in the com¬ 
plexity of the theory. As such, a different path will be 
taken here. 

To go beyond the RPA for the multi-component model, 
Eq.(II.8), the range optimization scheme will be used. 
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This scheme has been described in detail elsewhere, dl 
llfij but a summary is given here. 

The aim is to approximate the higher order terms em¬ 
bodied in Uij(k,Lu) in some manner. First, let Uij(k,ui) 
be independent of frequency to. This approximation is 
not necessary, but is a sensible one for computing the 
static equilibrium properties of the gas, such as gij{r) 
and E corr . Next, let Uij{k ) be real. This follows the com¬ 
mon assumption that the static local field factor G+(k) 
is dominated by its real part. [6] In this manner, the in¬ 
verse Fourier transform of Uij(k), Uij(r), can be viewed 
as a short-ranged attractive potential that counteracts 
the strong electron-electron Coulomb repulsion in v(r) 
at small r. What then is Uij(r )? 

Now, the RPA works well at high density, r s < 1, 
where the kinetic energy and exchange interactions dom¬ 
inate the Coulomb repulsion. At low density, r s 
1, where the RPA breaks down, the electron-electron 
Coulomb repulsion causes gij(r ) to be essentially zero 
out to some range cr,; 7 (see Figure 1). Notice though 
that if gij (r) were zero, it makes little difference to the 
electrons that the repulsion at that distance were infi¬ 
nite as opposed to just very large. As such, replace the 
Coulomb potential with a hard-core one for distances 
r < <Jij. While the RPA fails for hard-core potentials nat¬ 
urally, methods developed in the classical theory of liq¬ 
uids have found ways to overcome this problem. |22| One 
is the mean spherical approximation (MSA) closure. [ 23 ] 

Applied to jellium, the MSA states that if a^ is the 
range of the hard-core potential between electrons of spin 
i and j, then gij(r) is zero inside this range. But since 
gij (r) is zero inside, the equations relating (r) to gij (r) 
can be used to determine Uij{r ) inside. That is, Uij(r) 
takes whatever form is needed to ensure that gij{r ) is 
zero inside the core. The closure is summarized as: 


gij{r)= 0, r < a tj , 
Uij(r)= 0 , r > 


(II. 10) 


This closure, along with Eqs. (II.5), (II.6) and (II.81, 
form a closed set of equations for gij{r) and Uij{r), as- 
suming the aij are known. The last step is to optimize 
the range, aij, by letting it have the smallest value such 
that g^ (r) is positive for all r. Since the theory will now 
work properly for low and high densities, it presumably 
will work well for intermediate densities as well. This set 
of self-consistent equations will be referred to as RO-RPA 
theory. 

There are at least two ways to compute the compress¬ 
ibility: by the structure route using y(fc, 0) in Eq. (II.3), 


and by the energy route using an expression for the total 
energy[Q. Since the RO-RPA theory is approximate, the 
structure and energy routes will not give the same value 
for K. It is well known though that enforcing consis¬ 
tency between these two routes, i.e., enforcing the “sum 
rule”, can improve a theory, often greatly. [24] This ther¬ 
modynamic consistency can be attained under range op¬ 
timization by noticing that g t] (r) need not be set to zero 
inside the core. Instead, the minimum of gij(r) could 


have a non-zero value, go. The MSA closure can then be 
generalized to: 

9 ij{r)= 90, r < aij , 

Uij( r )= 0, r > a^, ( 11 . 11 ) 

where the value of go is determined by enforcing the sum 
rule on the compressibility. This extra condition, along 
with the RO-RPA equations above, will be referred to as 
“thermodynamically consistent” RO-RPA (TCRO-RPA) 
theory. 


III. NUMERICAL SOLUTION 


All theories were solved numerically as follows. Func¬ 
tions of r or k were solved on a grid of N r points 
with spacing Ar or Afc = ir/(N r Ar), respectively. Un¬ 
less stated otherwise, N r and Ar were set to 2 11 and 
0.05 /kF, respectively. To compute the static structure 
factor Sg(k ), the integration of Xij(k,ui) over to given 
by Eq. (II.5) was performed along the positive real axis. 
It is well known that along the real axis, a contribu¬ 
tion to Sij(k) from the plasmon mode must be accounted 
for 0(25] and that was done. As a check though, the inte¬ 
gration was also performed along the positive imaginary 
axis. 01115 ] In either case, the integral over frequency was 
evaluated using Romberg integration[5Bj with the rela¬ 
tive error tolerance being 10 -6 and 10~ 9 for the real and 
imaginary axis cases, respectively. 

O nce Sg(k) was computed, the integral over k in 
Eq.( |U.6 ) was evaluated by inverse Fourier transform to 
obtain g,j (r). Asa check on the accuracy, the RPA values 
for gij{ 0) were computed. It was found that unlike the 
correlation energy, E corr , g t j( 0) was sensitive to the grid 
spacing, but setting N r = 2 16 and Ar = 0.05/(32fcp) 
gave convergent RPA values of gij( 0) within 0.1%.[27] 
The grid spacing did not affect greatly any other quan¬ 
tity, though for increased accuracy in determining the 
density at which the spin susceptibility \S diverged, N r 
and Ar were set to 2 13 and 0.0125 //cf, respectively, for 


r s > 70. 

A new algorithm was used to implement range opti¬ 
mization. This new algorithm has two advantages over 
the one used in past work jl5IfT8] : it is straightforward 
to apply to inhomogeneous liquids, and is more efficient 
even for jellium. It is as follows: 1) An initial guess is 
made for the optimized potentials Uij(r), which could 
be zero. 2) With the Uij(r), the theory is solved as de¬ 
scribed above for the pair correlation functions gij(r). 3) 
The difference A gij{r) = gij{r) — gO is computed for all 
r, where gO = 0 for standard range optimization. 4) As a 
variation on Picard iteration [22], the change in the value 
of the optimized potential is set to a A gij(r), with the 
mixing parameter a < 0.25. 5) Since the optimized po¬ 
tential is an attractive potential, its new value for each 
r is checked to determine if it is greater than zero; if so, 
it is set to zero. 6) The difference between the new and 
old values is checked to determine if the potential has 
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converged; if not, the steps starting at 2) are repeated 
until convergence is obtained. Here, the relative error 
tolerance on each point of Uij(r) was 10 -4 , although in 
some cases it was reduced to 10” 5 as a check. 

Note that in this algorithm the optimized ranges, <7,j, 
are not considered explicitly. Instead, the algorithm re¬ 
lies on knowing only the value of the pair correlation 
function, gij , to obtain a refined guess for the optimized 
potential, Uij, at the same point. In that way, the above 
algorithm can be used for inhomogeneous liquids by the 
mere replacement of gij(r) and u^ (r) with g^ iq, r 2 ) and 
u t .j (rj, r 2 ), respecti vely, and then using the inhomoge¬ 
neous analog of Eq.( |IL8| ).[25] 

For all theories except TCRO-RPA the correlation en¬ 
ergy was computed in the usual manner as a charg¬ 
ing integral over the coupling constant e 2 , i.e., via the 
adiabatic-connection fluctuation-dissipation theorem. [T9] 
Define e c = 2aoE corr /(Ne 2 ) as a scaled correlation en¬ 
ergy per particle (in units of Rydbergs), with N being 
total number of electrons. Then this energy equation is: 

e c {r s ) = — f dX w( A), (III.l) 

ar s J o 

where 

w(X) = -dk[S(k,X)-S(k,0)]. (III.2) 


Here, S(k, A) = 1/2 JV. Sij(k, X) is the total structure 
factor for a gas with electron-electron potential Xv(r). 
The integral over At, Eq.(III.2), was computed via quadra¬ 
ture. The charging integral, Eq.(III.l I, was computed us¬ 
ing Romberg integration with an error tolerance of 1CU 4 . 

Enforcing thermodynamic consistency for the TCRO- 
RPA theory was done as follows. First, for a given value 
of r s , a value for gO was guessed. Then the optimized 
potentials, rtij(r), were computed in the same manner as 
for the RO-RPA. Next, the compressibility was computed 
using the structure route formula, Eq. (II.31. For the com¬ 
pressibility via the energy route, its value was computed 
initially using the Perdew-Wang fit for the correlation 
energy E corr ,^Z§\ and an expression relating the total en¬ 
ergy to the compressibility. [B] The change in the value of 
g 0 was set proportional to the difference between values 
of K 0 /K from these two routes. With this new gO , steps 
2)-6) above were repeated, and this iteration was contin¬ 
ued until the value of gO converged. This procedure was 
repeated to obtain Xij(k,ui) on a grid for 0 < r s < 11 
(see below) with spacing A r s = 0.1. 

These density-density response functions were then 
used to compute the correlation energy over this range. 
The representation of E corr expressed as an integral over 
r s was used: HUE] 


e c (r s ) = e r c pa ( r s) - dxAj(x), (III.3) 

irotrj J Q 

where 

A 7 (r s ) = -xj— [ dk[S(k,r s ) - S rpa (k,r s )] . (III.4) 
J o 


Here, e rpa (r s ) and S rpa {k,r s ) are the RPA values for 
the correlation energy and total structure factor, respec¬ 
tively, at mean separation r s . An accurate interpolation 
formula for e rpa due to Perdew-Wang[25] was used here. 

The integral, Eq. (III.4), for Ay(r s ) w as com puted for 
each r s in the same manner as for Eq.(III.2). The en¬ 
ergy route expression for the compressibility consists of 
derivatives of e c with respect to r s . To minimize errors 
then, this set of Ay values was then fit to an nth degree 
polynomial in r s . Degree n = 11 was found to give a 
sufficiently accurate fit (n = 9 worked almost as well). 
With this functional form, Eq.(III.3) was evaluated ana¬ 
lytically. The self-consistent theory for g 0 was then solved 
again, with the new and old values for e c being used to 
compute the energy route K with a mixture of 1:1 old 
to new. After new density-density response functions, 
Xij(k,oj) were computed, the procedure to compute new 
values for e c was repeated. It was found that the fitted 
values for K had converged to within 10~ 3 (10~ 5 for e c ) 
for r s < 10 after seven iterations. 

At r s « 0, go ~ 0 naturally, then go rose to a maxi¬ 
mum of 0.177 at r s ss 1.7, and then gradually dropped to 
zero again at r s « 10.8. When the positivity constraint 
on gij(r ) was relaxed, go became slightly negative as r s 
increased beyond 10.8. Since this positivity constraint 
on gij (r) is more important than enforcing a sum rule, 
r s ~ 10.8 is then the limit of the usefulness of enforcing 
thermodynamic consistency on the compressibility. How¬ 
ever, it will be shown below that since the RO-RPA is 
most accurate at low density, this limit is not regarded 
as important. 

For comparison, some results of the theories of Singwi- 
Tosi-Land-Sjolander (STLS), [SB] and Utsumi and Ichi- 
maru (UI) [5fT will also be shown. The UI theory is con¬ 
sidered accurate for the short range behavior of the pair 
correlation function at metallic densities. [7j The STLS 
theory is considered accurate for the correlation energy 
and almost as accurate as UI for the pair correlation 
function, but is also straightforward to implement. The 
theory has also been generalized to apply to inhomo¬ 
geneous liquids in atoms and ions. [32] The STLS the¬ 
ory was solved for the total structure factor S(k) in a 
similar manner to that of Sij(k) for the RO-RPA. Once 
S(k) was determined self-consistently, the spin-averaged 
g(r) = 1/4 "}2ij9ij( r ) was obtained by inverse Fourier 
transform using the analog of Eq. (II.6). The UI theory 
was solved in the same manner as the RPA, but with 
v{k) v(k)( 1 — G + (k)). Values for the static local field 
factor G+(k) were interpolated from data presented in 
Table I of [3]]. 


IV. RESULTS 

Unless noted otherwise, all RO-RPA results given here 
will be for the multi-component version described in 
Sec|TT] above. Results of the one- and multi-component 
versions of the RPA are the same for the quantities ex- 
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amined here. 


Figure [T] shows results for the spin-averaged pair corre¬ 
lation function, g(r), for various densities. As a compar¬ 
ison, QMC simulation data of Ceperley and co-workers 
is also shown. [33l 34, As can be seen, the predictions of 
the RO-RPA are much improved over the RPA, with the 
RO-RPA outperforming, as expected, the STLS theory 
at very low density, r s = 50. For r s = 2, the contact 
value, <?(0) = 0.176 and 0.175 for the UI and TCRO- 
RPA theories, respectively, which are 10% smaller than 
the simulation value. For r s = 10, the thermodynam¬ 
ically consistent value of go was close to zero. Conse¬ 
quently, the TCRO-RPA prediction for g(r ) (not shown) 
is almost the same as the RO-RPA for this density (and 
lower densities). Holm and von Barth have shown that 
the one-shot and fully self-consistent GWA produce only 
modest improvement in the local structure of g(r) over 
the RPA, for the metallic density r s = 4. [531 

As mentioned above, the focus on static properties 
partly justified ignoring the frequency dependence of the 
optimized potentials, Uij. However, the structure of the 
multi-component theory is such that when mapped to 
the one-component form for x(fc,w), Eq.(II.9l, the local 
field factor G + that arises is frequency dependent. It is 
interesting then to examine theoretical predictions for a 
dynamic property of the gas: its collective excitations, 
i.e., plasmons. mm 


Figure [2] shows the plasmon frequency lo p as a function 
of wavevector k for various theories and simulation data 
for the same densities as given in Figure [I] As shown, 
the plasmon frequency is normalized by its value at zero 
wavenumber: w p (0) = 47me 2 / m. The curves termi¬ 
nate at the beginning of the electron-hole continuum of 
the non-interacting electron gas. [6] The Corradini et al. 
curves were obtained using their fit® [37] of G + (k, 0) to 
QMC simulation data of Moroni et al.[3H f° r r s = 2, 5 
and 10. While not shown in the figure, the predictions of 
the theory of UI[31] for r s = 2 are essentially the same 
as those shown for Corradini et al. As can be seen, the 
RO-RPA predicts larger plasmon frequencies than the 
STLS theory for all densities, with the difference increas¬ 
ing somewhat as r s increases. The RO-RPA predictions 
agree well with the results using the Corradini et al. fit 
for both densities, r s = 2 and 10. It has been shown 
elsewhere that the fully self-consistent GWA gives poor 
predictions for the spectral properties of jellium, includ¬ 
ing the plasmon modes. [H] 


Figure [3] shows theoretical predictions for the scaled 
correlation energy per particle, e c , as a function of the 
scaled average electron separation, r s . Also shown are 
simulation data of Ortiz et al. [32], and the Perdew-Wang 
fit [29] to simulation data of Ceperley and Alder [33], 
Vosko et al.[20] show results for e c for other theories for 
the paramagnetic phase. As can be seen, the predictions 
of the RO-RPA greatly improve upon those of the RPA. 
As mentioned above, the density range of most metals is 
2 < r s < 6. [4D] For this range, the RO-RPA values for 
e c are more negative than those of simulation by 15% on 
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FIG. 1. Spin-averaged pair correlation function g(r) as a 
function of distance r for various densities. The meanings of 
the curves and symbols are shown in the figure legend. For 
r s = 2 and 10, the simulation data are from Ceperley and 
Alder[331 as obtained from Gori-Giorgi et al.[36] For r 3 = 50, 
the simulation data are from Zong et al. [34] The contact value 
g( 0) from the RPA theory is -0.66, -3.95 and -15.1 for r s = 2, 
10 and 50, respectively. 




average, while the RPA values are more negative by 46% 
on average. The accuracy of the RO-RPA increases with 
decreasing density such that, for example, at r s = 40, its 
value for e c is within 3% of simulation. The thermody¬ 
namically consistent theory, TCRO-RPA, improves even 
more on the RPA: its predictions for e c are within 4% 
of simulation over the metallic density range. Interest¬ 
ingly, the predictions of the one-component RO-RPA the¬ 
ory, while being inferior to the multi-component theory 
for almost every other quantity, are slightly better than 
the multi-component for e c when thermodynamic consis¬ 
tency is enforced. The one-component TCRO-RPA the¬ 
ory agrees with the Perdew-Wang fit to within 3% for the 
metallic density range. This accuracy is the same or only 















6 



k/k F 

FIG. 2. Plasmon dispersion curves for various scaled average 
electron separations r s . The meanings of the curves are shown 
in the figure legend. An r s label refers to the closest RO-RPA 
curve. For the other theories, the trend is for the slope of 
cjp(k) at small k to decrease as r 3 increases. The complete 
RPA curve for r s = 50 is not shown, it terminating with 
value Ld p (k)/tjjp(0) = 1.55 at k/kF = 2.36. The Corradini 
et al. curves were obtained using their fit [37] of G+(fc,0) to 
simulation data[55j. 



FIG. 3. Scaled correlation energy per particle e c (in Ryd¬ 
bergs) as a function of the scaled average electron separation 
r s . The meanings of the curves and symbols are shown in 
the figure legend. The Perdew-Wang curve is a fit 125] to 
simulation data 133) . 



FIG. 4. Scaled inverse compressibility Ko/K as a function 
of the scaled average electron separation r s , obtained using 
Eq.( |II.3[ ). The meanings of the curves are shown in the figure 
legend. Note that the TCRO-RPA and Perdew-Wang curves 
overlap for almost all densities shown. 



FIG. 5. Scaled inverse spin susceptibility Xp/xs as a func¬ 
tion of the scaled average electron separation r s , obtained 
using Eq.( |IL4| ). The meanings of the curves are shown in the 
figure legend. 


slightly less than that of the most accurate theories for 
e c : STLS,® and Vashishta and Singwi.flTI which agree 
with simulation within 2% and 3%, respectively, over 
this density range (the UI theory agrees within 7%[5T]). 
While it has been tested to date for only a few densities, 
the fully self-consistent GWA gives very good agreement 
for e c , within 1% for r s = 2 and 4- nn Given the mediocre 
predictions of the theory for other properties, this good 
agreement is thought to be to due to a large cancellation 
of effects. [551 

Figure [4] shows results for the scaled inverse of the 
compressi bilit y, Kq / K , as defined by the structure equa¬ 
tion, Eq. W ) above. Simulation values were obtained 
using the Perdew-Wang fit [25] to data of Ceperley and 
Alder [33 for the correlation energy, and an expression 
relating the compressibility to the total energy[6]. For 
jellium, simulation studies[55] show that the inverse com¬ 
pressibility goes to zero at r s = 5.25. The RO-RPA pre- 
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diets a zero at r s = 6, which is 14% higher. As a compar¬ 
ison, the STLS theory predicts a zero at r s = 3, which 
is 55% less than simulation. As is well known, the RPA 
gives the non-interacting value at all densities. Interest¬ 
ingly, the TCRO-RPA predictions agree very well with 
the Perdew-Wang fitted data, within 0.1%, up to r s = 8. 


Figure [5] shows results for the scaled inverse of the 
spin susceptib ility, \p/Xs > as defined by the structure 
equation, Eq.(II.4| above. As for K , simulation val¬ 
ues were obtained using the Perdew-Wang fit [25] for 
the correlation energy, and an expression relating xs 
to the total energy [ 5 ] ■ A divergence in \S is identi¬ 
fied with a second-order transition from the paramag¬ 
netic state to a polarized state. In Hartree-Fock theory, 
this polarized state was identified with the fully polar¬ 
ized, i.e., ferromagnetic, state. [19] Simulation work ap¬ 
peared to have reinforced this idea. [33] In the Hartree- 
Fock and RPA theories, the density of the divergence 
of xs is slightly lower (higher r s ) than the density for 
paramagnetic-ferromagnetic phase coexistence. In their 
fit then, Perdew and Wang reasoned that the true den¬ 
sity of divergence occurs at an r s slightly above that for 
paramagnetic-ferromagnetic coexistence, r s = 73, which 
was computed via simulation by Ceperley and Alder. (33 
This reasoning yielded a divergence of xs at r s = 77.5. 
This zero of 1 /xs can be seen in Figure [5j 

Subsequent simulation work by Ortiz et al. found 
that partially polarized states were energetically favor¬ 
able at higher densities (lower r s ), and the transition 
from the paramagnetic to a partially polarized state was 
continuous. |39| Later simulation work by Zong et al. un¬ 
derscored this picture, though their estimate of the den¬ 
sity of this transition was lower, r s ss 50. jM] In addition, 
Zong et al. estimated the transition to the ferromagnetic 
state was at r s « 100, which was also their estimate for 
the transition to the Wigner solid state. 


As can be seen in Figure [5j the RO-RPA predicts a 
divergence at r s ~ 107, in good agreement with the esti¬ 
mate of Zong et al. for the paramagnetic-ferromagnetic 
transition. So, either the RO-RPA overestimates the 
value of r s for the paramagnetic-partially polarized tran¬ 
sition by a factor of two, or predicts that only the 
paramagnetic-ferromagnetic transition is second order. 
Given the accuracy of the RO-RPA for the correlation en¬ 
ergy at low density, the latter explanation appears more 
plausible. However, the RO-RPA value for the density 
of divergence was obtained from the structure route ex¬ 
pression for the susceptibility, and usually estimates us¬ 
ing the energy route are more accurate. Given the from 
of the RO-RPA multi-component theory, it is straightfor¬ 
ward to generalize it to compute E corr as a function of 
the fractional spin polarization, p = (n-f(x) — n±(x))/n, 
and thus determine coexistence and spinodal boundaries. 
The solution to this puzzle then is left to future research. 

In this work, thermodynamic consistency was enforced 
on the compressibility, which is a measure of the sensi¬ 
tivity of the electron density to changes in the pressure. 
It is not expected then that this method would neces¬ 


sarily improve the predictions of the theory for the spin 
susceptibility, which is a measure of the sensitivity of a 
different quantity, the spin polarization, to changes in 
the magnetic field. Nonetheless, it is interesting to ex¬ 
amine the TCRO-RPA predictions for xs■ As can be 
seen in Figure [5] the TCRO-RPA values agree well with 
the Perdew-Wang fit at very high density (small r s ), but 
drop below the Perdew-Wang curve at r s ~ 0.5. For 
example, at r s = 2, the TCRO-RPA value is 10% below 
the Perdew-Wang one. For r s > 2, the TCRO-RPA curve 
asymptotically approaches the RO-RPA curve, terminat¬ 
ing at it at r s « 10.8. At that point, the TCRO-RPA 
curve is above the Perdew-Wang curve. In other words, 
the agreement with the fitted simulation data appears 
to be better for the RO-RPA than for the TCRO-RPA. 
It has been remarked elsewhere that the Perdew-Wang 
fit is not that accurate for polarized states. [6] It seems 
unlikely though that this inaccuracy would be that large 
near p = 0. A clarifying task then would be to com¬ 
pute xs via the energy route. For this, E corr would be 
computed for small p and constant r s , while enforcing 
thermodynamic consistency on K (or even xs )■ 


V. SUMMARY AND DISCUSSION 


In summary, the range optimization scheme was ap¬ 
plied to the RPA theory for jellium. It was shown that 
this RO-RPA theory gives greatly improved predictions 
for the gas properties as shown by its results for the pair 
correlation function g(r), compressibility K and spin sus¬ 
ceptibility xs- For the correlation energy, E corr , the the¬ 
ory is most accurate at low densities, but it still gives pre¬ 
dictions within 15% of simulation for the density range of 
most metals, 2 < r s < 6. Enforcing thermodynamic con¬ 
sistency on the compressibility improves the agreement 
with simulation to within 4% for this density range, while 
the one-component version of the theory is slightly better 
at 3%. This agreement is comparable to the most accu¬ 
rate of the previous theories for E corr .[ 14, !30, 144] Also, 
the RO-RPA appears to outperform the STLS theory in 
comparison with simulation data for the plasmon modes 
and compressibility, and the pair correlation function at 
low density. 

The thermodynamically consistent theory could be fur¬ 
ther improved by conducting the range optimization to 
obey the cusp condition on g(r) near r = 0.02] Also, 
since range optimization can be applied to any theory 
in which the positivity condition of g(r ) is violated, im¬ 
provements in the accuracy of the basic theory are pos¬ 
sible. 

One noteworthy result was that the RO-RPA theory 
predicts a divergence of xs at a density that is lower than 
current estimates for the liquid-solid transition. 0 [M] 
Thus, no divergence of xs is predicted for the liquid 
phase. This result was obt ained using the structure 
route expression for xs , Eq.(II.4|. Given the evidence 


from simulations for a transition to partially polarized 




states, [HJ [35] it would be interesting within the RO-RPA 
to examine the phase behavior of jcllium as a function of 
polarization using the energy route. m Further, since the 
temperature dependent density-density response func¬ 
tion can be represented by an equation almost identical 
to Eq.(II.7),pj)j the full temperature dependent phase di¬ 
agram can be obtained. 

As stated above, one aim of this work is to apply range 
optimization to inhomogeneous electron liquids, to com¬ 
pute the band structure of semiconductors, for example. 


It was shown in Sec III how the new algorithm to im¬ 
plement range optimization can be applied to inhomoge¬ 
neous liquids, and results for that will be presented in a 
future work. 
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